{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "739bad01",
   "metadata": {},
   "source": [
    "## 1.2.2 常见离散型分布\n",
    "\n",
    "#### 二项分布"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d59b8794",
   "metadata": {},
   "source": [
    "> 例子：\n",
    "> 现有90台同类型的设备，各台设备的工作是相互独立的，发生故障的概率是0.01，且一台设备的故障能由一人处理，配备维修工人的方法有两种，一种是3人分开维护，每人负责30台，另一种是由3人共同维护90台，试比较两张方法在设备发生故障时不能及时维修的概率的大小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f2147b06",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "58fedc80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. 3人分别负责30台，发生故障时不能及时维修表明30台机器中发生故障的机器台数为2台及以上的概率，由于有3组，在计算出魅族的概率后需乘以3\n",
    "# 2. 3人共同维护90台机器，则是90台机器中同时发生故障的机器超过4台及以上的概率是多少\n",
    "# 3. 可以使用二项分布来求概率，由于计算的是随机变量大于某个值的概率分布，因此需使用生存函数进行计算（《生存分析与可靠性》）\n",
    "# cdf + sf = 1\n",
    "# cdf：累加分布函数\n",
    "# sf：生存函数\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "010de1b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p30 =  0.10844399493672194 , p90 =  0.012940166077704174\n",
      "3人共同维护90台机器时，发生故障不能及时修理的概率更小。\n"
     ]
    }
   ],
   "source": [
    "# 通过生存函数计算发生故障不能及时维修的概率\n",
    "\n",
    "# 3人分别负责30台发生故障不能维修的概率\n",
    "# k 成功次数, p 发生故障的概率, n 实验的总次数\n",
    "p30 = st.binom.sf(k=1, p=0.01, n=30) * 3\n",
    "\n",
    "# 3个人共同负责90台设备发生故障不能修理的概率\n",
    "p90 = st.binom.sf(k=3, p=0.01, n=90)\n",
    "\n",
    "print('p30 = ', p30, ', p90 = ', p90)\n",
    "print('3人共同维护90台机器时，发生故障不能及时修理的概率更小。')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3eba1e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 通过概率质量函数和累计分布函数计算概率\n",
    "\n",
    "# 还可以通过另一种方式计算3组每组30台机器发生故障不能及时修理的概率：\n",
    "# 1. 计算30台机器不发生故障的概率与同时只有1台发生故障的概率之和。\n",
    "# 2. 然后再用1减去该概率和，即是同时发生2台及以上故障的概率\n",
    "# 3. 最后还要乘以3，表示共3组每组30台机器发生故障不能及时修理的概率\n",
    "# 4. 3人共同负责90台发生故障不能及时修理的概率的计算类似\n",
    "# 5. 因此离散型随机变量的生存函数概率计算等于P(x>k)，即不包含k值本身，而累积分布函数的概率计算等于P(x<=k)，是包含了k值的\n",
    "# 6. 离散型随机变量的概率质量函数P(x=k)是有意义的，连续型随机变量的概率密度函数值是无意义的，P(x=k) = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4033f278",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.10844399493672308\n",
      "0.10844399493672208\n"
     ]
    }
   ],
   "source": [
    "# 下面两种计算方式的效果相同，pmf函数计算离散型随机变量的概率质量函数\n",
    "p30_1 = (1-(st.binom.pmf(k=0, p=0.01, n=30)+st.binom.pmf(k=1, p=0.01, n=30))) * 3\n",
    "\n",
    "# cdf 计算累积分布函数的值\n",
    "p30_2 = (1-st.binom.cdf(k=1, p=0.01, n=30)) * 3\n",
    "\n",
    "print(p30_1)\n",
    "print(p30_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "596a0afe",
   "metadata": {},
   "source": [
    "#### 泊松分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0ea6727b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3000英尺的金属线上有5个及以上疤痕的概率是（通过累积分布函数）： 0.185\n",
      "\n",
      "3000英尺的金属线上有5个及以上疤痕的概率是（通过生存函数）： 0.185\n"
     ]
    }
   ],
   "source": [
    "# 1. P(x>=5)的概率等于k=5到无穷的概率质量之和。可以用1-P(x<=4)来间接计算。\n",
    "#    这里用到Scipy统计包的cdf函数，即“累积分布函数（Cumulative distribution function）”\n",
    "# 2. 直接通过生存函数sf 也可以求 P(x>=5)。sf是“生存函数（survival function）”\n",
    "# 3. cdf 和 sf 两个函数后面会经常用到。对于相同分布的相同随机变量，cdf + sf = 1\n",
    "\n",
    "# 通过累积分布函数间接求解\n",
    "# mu 表示均值，本例中等于3，下同。\n",
    "p4 = st.poisson.cdf(k=4, mu=3)\n",
    "p5 = 1 - p4\n",
    "print('3000英尺的金属线上有5个及以上疤痕的概率是（通过累积分布函数）：', np.round(p5, 3))\n",
    "\n",
    "# 通过生存函数直接求解\n",
    "# 注意这里的k不是5，而是4\n",
    "p5_1 = st.poisson.sf(k=4, mu=3)\n",
    "print('\\n3000英尺的金属线上有5个及以上疤痕的概率是（通过生存函数）：', np.round(p5_1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fece54dc",
   "metadata": {},
   "source": [
    "## 1.2.3 常见连续型分布\n",
    "\n",
    "#### 正太分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a989b185",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "正态分布概率密度函数曲线（均值：1，标准差：0.5）：\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnaklEQVR4nO3deXRU93338fd3RrvQBhKbNjYZI5vNyIB3J16CbWrsrLhJ46Rx/dDGaXrapHFO25ynJ33S5EmXtI1Th+O6SU8WP06cJjgmdpzENrjGgEACA2IRAi0IkNC+SzPzff6YEZaFhEbSjO4s39c5Osyd+5s737Glj65+93d/P1FVjDHGRD+X0wUYY4wJDQt0Y4yJERboxhgTIyzQjTEmRligG2NMjEhw6o1zc3N10aJFTr29McZEpQMHDlxS1byx9jkW6IsWLaK8vNyptzfGmKgkIrXj7bMuF2OMiREW6MYYEyMs0I0xJkZYoBtjTIywQDfGmBgRVKCLyCYROSEi1SLy5Bj7s0TkRRE5JCJHReTToS/VGGPM1UwY6CLiBp4C7gNKgUdEpHRUs88Cx1R1NXAn8I8ikhTiWo0xxlxFMOPQ1wPVqloDICLPAVuAYyPaKJAhIgLMAloBT4hrNSasVJW9Z1o5VN9OapKb20vyWJSb7nRZxgQtmEDPB+pHbDcAG0a1+TawA2gEMoCPqapv9IFE5HHgcYCioqKp1GtMWNS39vL55yo4WNd++TkR+PiGIv76gVJSEt3OFWdMkIIJdBnjudGrYnwAqATeDywFXhWR3ara+Z4XqW4HtgOUlZXZyhomIlQ3dfOx7+5h0Ovj7z+4kvtXLqCrf4hndp/he2+d5XRTD//56Rst1E3EC+aiaANQOGK7AP+Z+EifBn6mftXAGeDa0JRoTPi09w7y6LP7EIGff/YWHllfRFZqIgU5afzvB6/jnz66mrfPtPDFnx7GVvcykS6YQN8PlIjI4sCFzq34u1dGqgPuAhCRecByoCaUhRoTaqrKky+8Q1NXP//x6I0szZt1RZsP3lDAF+5dzouHGvnJgQYHqjQmeBMGuqp6gCeAV4Aq4HlVPSoi20RkW6DZV4GbReQd4LfAl1T1UriKNiYUXjvRxMtHL/Dn9yxndWH2uO3++I6lbFg8m7/75TFaugdmrkBjJkmc+jOyrKxMbbZF45RBj49N39oFwMt/djtJCVc/t6lu6mLTt3bzkbJC/v6DK2eiRGPGJCIHVLVsrH12p6iJSz/eV0fNpR7+ZnPphGEOsGxuBp/YWMxPyuupb+2dgQqNmTwLdBN3hrw+tu+qYf2i2bzv2rlBv27bHUtxifCd10+HsTpjps4C3cSdXx5u5Fx7H9vuXDKp183PSuFjNxby0wP1NLb3hak6Y6bOAt3EFVXlu2/UsHxeBu9bHvzZ+bD/dccSvD7lR3vrwlCdMdNjgW7iSnltG8cvdPGZWxfjn6licgpy0nj/tXN5bn89g54rboY2xlEW6CauPLevnlnJCWxevWDKx/j4xmIudQ/w62MXQliZMdNngW7iRmf/EC+908iDaxaSljT19dHvKMmjcHaqdbuYiGOBbuLGjspG+od8bL2xcOLGV+FyCR+6oYA9NS1c6OgPUXXGTJ8Fuokbv6g8R8ncWazMz5r2sR5ak48q7Dh0LgSVGRMaFugmLpzv6GP/2TYeXL1wShdDR1uUm87qwmx+XjF6njpjnGOBbuLCS4fPA7B59cKQHfOhNQs5dr6TUxe7QnZMY6bDAt3EhRcPn+e6hZksDuEKRA+sWoAIvHzERruYyGCBbmLeufY+DtW3s3lV6M7OAeZmpLCmMJtXqy6G9LjGTJUFuol5vw0E7r3XzQv5se8tnc/hhg6bCsBEBAt0E/N+U9XE4tz0MRewmK57SucF3sPO0o3zggp0EdkkIidEpFpEnhxj/xdFpDLwdUREvCIyO/TlGjM53QMe3j7dwl2TmFVxMpbNncWSvHRePWaBbpw3YaCLiBt4CrgPKAUeEZHSkW1U9ZuqukZV1wBfBt5Q1dYw1GvMpLx5qplBr4+7VoS+u2XYPaXz2HO6hY6+obC9hzHBCOYMfT1Qrao1qjoIPAdsuUr7R4Afh6I4Y6brN1VNZKYkULYoJ2zvcW/pPDw+5Y2TzWF7D2OCEUyg5wP1I7YbAs9dQUTSgE3AC+Psf1xEykWkvLnZvvlNeHl9ymvHm7hz+VwS3eG7XLSmMIfMlAR2W6AbhwXzXT7WbXXjLUT6e8D/jNfdoqrbVbVMVcvy8vKCrdGYKTlyroOWnkHuWhGe/vNhbpdwa0kuu04149QavcZAcIHeAIyczagAGO9+561Yd4uJEG9WXwLglmW5YX+v20vyuNg5wKmm7rC/lzHjCSbQ9wMlIrJYRJLwh/aO0Y1EJAu4A/hFaEs0Zmp2n2pmxYJMcmclh/29brvG/xfnLut2MQ6aMNBV1QM8AbwCVAHPq+pREdkmIttGNH0Y+LWq9oSnVGOC1zvo4WBtO7eVhP/sHCA/O5Wleel2YdQ4KqhZ/lV1J7Bz1HNPj9r+HvC9UBVmzHTsO9PKoNfHrTPQ3TLs9mvy+NHeOvqHvKQkumfsfY0ZZneKmpj05qlLJCW4WL945u5vu70kjwGPj31n7BYM4wwLdBOT3qy+RFlxzoyeKW9YMpsEl/B2TcuMvacxI1mgm5jT3DXA8Qtd3DpD/efD0pISWFWQZYFuHGOBbmLO3jP+QL156cwGOsDGJXM43NBBz4Bnxt/bGAt0E3P2nWklLcnN9QszZ/y9NyyZg8enHKhtm/H3NsYC3cScvTWtrCvOISGMt/uPp6w4B7dLLv+VYMxMskA3MaWtZ5ATF7vYuGSOI++fnjzcj24jXczMs0A3MWXfWX+QzuRwxdE2LpnDofp2egetH93MLAt0E1P21rSSnOBiVUGWYzVstH504xALdBNT9p1tYW1RNskJzt2puW64H926XcwMs0A3MaOzf4hjjZ1sWOxM//mwWckJrMzPYo+NRzczzALdxIwDZ9vwKWxwsP982IbFs3mnoYP+Ia/TpZg4YoFuYsbbZ1pIdAtri8K33Fyw1hXnMOj1ceRch9OlmDhigW5ixoGzbazMzyI1yfmZDm8o9v9SKbcLo2YGWaCbmDDk9fHOuY6IODsHyJ2VzOLcdBvpYmaUBbqJCcfPdzHg8bG2KNvpUi67oSiHg7Vtts6omTFBBbqIbBKREyJSLSJPjtPmThGpFJGjIvJGaMs05uoq6v1nwmsKs50tZISyRTm09AxytqXX6VJMnJgw0EXEDTwF3AeUAo+ISOmoNtnAd4AHVfU64COhL9WY8VXWtZOXkUx+dqrTpVy2brgf/ayNRzczI5gz9PVAtarWqOog8BywZVSb3wd+pqp1AKraFNoyjbm6ivp21hRmIyJOl3LZsrxZZKYkcLDO+tHNzAgm0POB+hHbDYHnRroGyBGR10XkgIh8cqwDicjjIlIuIuXNzbaYrgmNtp5Bzlzqiaj+cwCXS7ihOIfysxboZmYEE+hjnfKMvsqTAKwDHgA+APyNiFxzxYtUt6tqmaqW5eXlTbpYY8ZS2dAOwNrCyBjhMlJZcQ6nmrrp6B1yuhQTB4IJ9AagcMR2AdA4RpuXVbVHVS8Bu4DVoSnRmKurrGvHJTg6Idd4hsejW7eLmQnBBPp+oEREFotIErAV2DGqzS+A20QkQUTSgA1AVWhLNWZsFfXtXDMvg/TkBKdLucKawmzcLrHx6GZGTPgToKoeEXkCeAVwA8+q6lER2RbY/7SqVonIy8BhwAc8o6pHwlm4MQA+n1JZ18YDqxY4XcqY0pISKF2QSXmtjXQx4RfUKY2q7gR2jnru6VHb3wS+GbrSjJnYmZYeOvs9Edl/PmxdcQ7/b389Hq/PkWXxTPyw7y4T1Srq2gFYE2EjXEZaW5RN35CXExe7nC7FxDgLdBPVKuvbyEhOYFneLKdLGdfw3auV9e2O1mFinwW6iWoVde2sKszC5YqcG4pGK5qdxuz0JCoDf00YEy4W6CZq9Q16OX6hK6L7zwFEhNUFWXaGbsLOAt1ErXfOdeD1acTdITqWtUU5VDd309VvNxiZ8LFAN1Groi7yZlgcz5rCbFThcIOtYGTCxwLdRK3K+naKZqcxZ1ay06VMaLVdGDUzwALdRK2KuvaoODsHyEpNZEle+uW/KowJBwt0E5XOd/RxobM/KvrPh60pzKayvt1WMDJhY4FuotLwEMBoOUMHWFuYzaXuQRra+pwuxcQoC3QTlSrr20lyuyhdmOl0KUFbExheaf3oJlws0E1Uqqhrp3RhJskJbqdLCdq1CzJITnBZoJuwsUA3Ucfj9XH4XHtU9Z8DJLpdrMy3G4xM+Figm6hz/EIX/UM+1hZF9h2iY1lTmM2Rcx0MeX1Ol2JikAW6iToVgTPctVF0QXTYmqJsBjw+jp+3mRdN6Fmgm6hTWddO7qwkCnJSnS5l0t6dedHGo5vQCyrQRWSTiJwQkWoReXKM/XeKSIeIVAa+vhL6Uo3xq6hvY01hNiKRO8PiePKzU8mdlXz5rwxjQmnCFYtExA08BdyDfzHo/SKyQ1WPjWq6W1U3h6FGYy7r6B2iprmHD91Q4HQpUyIi/huMbCpdEwbBnKGvB6pVtUZVB4HngC3hLcuYsVU2tAPRdUPRaGuLsqm51ENHr828aEIrmEDPB+pHbDcEnhvtJhE5JCK/EpHrxjqQiDwuIuUiUt7c3DyFck28q6xrRwRWFWQ5XcqUXe5HD/xyMiZUggn0sToqR09GcRAoVtXVwL8BPx/rQKq6XVXLVLUsLy9vUoUaA/7+82vmZpCRkuh0KVO2qiALEazbxYRcMIHeABSO2C4AGkc2UNVOVe0OPN4JJIpIbsiqNAZQVSrro2eGxfFkpCRSMneWjXQxIRdMoO8HSkRksYgkAVuBHSMbiMh8CQw5EJH1geO2hLpYE9/OtvTS3jsUdXeIjsVmXjThMGGgq6oHeAJ4BagCnlfVoyKyTUS2BZp9GDgiIoeAfwW2qn2nmhC7vEJRTAR6Dm29Q9S19jpdiokhEw5bhMvdKDtHPff0iMffBr4d2tKMea/K+nbSk9yUzM1wupRpG+42qqhrp3hOurPFmJhhd4qaqFFR186qgmzcrui7oWi0a+bNIjXRbRN1mZCyQDdRoX/IS9X5zpjoPwdIcLtYWZBld4yakLJAN1HhyLkOPD6N+hEuI60tyqaqsZMBj9fpUkyMsEA3UWG4ayIWLogOW1uYzaDXx7HGTqdLMTHCAt1EhYq6dgpyUpmbkeJ0KSFjS9KZULNAN1Ghoq4tprpbAOZnpTA/M4UKu2PUhIgFuol4Fzv7aezoj8oViiYyfIORMaFggW4i3vAZbKydoYP/mkBday8t3QNOl2JigAW6iXiV9e0kuoXrFmY6XUrIDf+SOmQzL5oQsEA3Ea+iro3SBZmkJLqdLiXkVhVk4bKZF02IWKCbiObx+njnXEdM9p8DpCUlsHx+pt1gZELCAt1EtJMXu+kd9MZk//mw4QujPp/NZ2emxwLdRLThESCxcsv/WNYWZtPV76HmUo/TpZgoZ4FuIlpFXRuz05Momp3mdClhM3z3qw1fNNNlgW4iWkVghaLA+ikxaWneLGYlJ9gKRmbaLNBNxOroG6K6qZu1Mdx/DuB2CasKsuwM3UxbUIEuIptE5ISIVIvIk1dpd6OIeEXkw6Er0cSr4YC7oTg2R7iMtLYom+Pnu+gbtJkXzdRNGOgi4gaeAu4DSoFHRKR0nHbfwL9UnTHTdrC2DRH/WO1Yt6YwB49POdLY4XQpJooFc4a+HqhW1RpVHQSeA7aM0e5zwAtAUwjrM3Gsor6d5fMyyEhJdLqUsBselmk3GJnpCCbQ84H6EdsNgecuE5F84GHgaa5CRB4XkXIRKW9ubp5srSaO+HxKRV1bzN5QNFpeRjL52anWj26mJZhAH2t4weg7IL4FfElVr9oBqKrbVbVMVcvy8vKCLNHEo9PN3XT1e2J6/Ploa4ps5kUzPcEEegNQOGK7AGgc1aYMeE5EzgIfBr4jIg+FokATn4ZnWLwhTs7QwX+D0bn2Ppq6+p0uxUSpYAJ9P1AiIotFJAnYCuwY2UBVF6vqIlVdBPwU+BNV/XmoizXx42BdG1mpiSzJTXe6lBkz/NeI9aObqZow0FXVAzyBf/RKFfC8qh4VkW0isi3cBZr4dDCwQpHLFbs3FI123cIsElxiE3WZKUsIppGq7gR2jnpuzAugqvqp6Zdl4lln/xCnmrp5YOVCp0uZUSmJblYsyLQzdDNldqeoiTiH6ttRje0JucazpjCbww3teG3mRTMFFugm4hysbUfk3Umr4snaomx6Br2cvNjldCkmClmgm4hTUd9GydxZZMbBDUWjlRXPBqC81ibqMpNngW4iiv+GonbWFsbPcMWRCmenkpeRzIGzrU6XYqKQBbqJKDWXeujoG+KG4mynS3GEiFBWnGNn6GZKLNBNRDlY5w+yeLnlfyxli2bT0NbHhQ67wchMjgW6iSgHzvpvKFqWN8vpUhxTFpguuLzWul3M5Figm4iyv7aVsuKcuLqhaLTShZmkJropP2vdLmZyLNBNxGjpHqCmuYeyRbOdLsVRiW4Xqwuz7AzdTJoFuokYwxcCb1wUv/3nw25cNJuq8130DHicLsVEEQt0EzH2n2klKcHFyjhYoWgi64pz8PrUptM1k2KBbiLG/to21hRkk5zgdroUx91QnIMI7Lfx6GYSLNBNROgd9HD0XAdl1t0CQGZKIsvnZXDAxqObSbBANxGhsr4dj0+5Mc4viI5UtiiHg7VteLw+p0sxUcIC3USE/WfaEPF3NRi/suLZ9Ax6OX7BJuoywQkq0EVkk4icEJFqEXlyjP1bROSwiFQGFoG+NfSlmlhWXtvK8nkZZKXG34Rc41kX+OVm3S4mWBMGuoi4gaeA+4BS4BERKR3V7LfAalVdA/wh8EyI6zQxzOP1cbC2zbpbRinISWV+ZopdGDVBC+YMfT1Qrao1qjoIPAdsGdlAVbtVdXhG/nTAZuc3QTt+oYueQS83LrZAH0lE2LBkNm/XtPLuj5cx4wsm0POB+hHbDYHn3kNEHhaR48BL+M/SjQnKvjP+M9Ay6z+/wsYlc7jUPcDp5h6nSzFRIJhAH2tSjStOF1T1v1X1WuAh4KtjHkjk8UAfe3lzc/OkCjWxa09NC0Wz01iYnep0KRHnpiVzAP9/I2MmEkygNwCFI7YLgMbxGqvqLmCpiOSOsW+7qpapalleXt6kizWxx+tT9ta0XA4u817Fc9KYn5nC2xboJgjBBPp+oEREFotIErAV2DGygYgsExEJPL4BSALsO9BMqOp8J539Hm5aaoE+FhHhpqVz2FvTYv3oZkITBrqqeoAngFeAKuB5VT0qIttEZFug2YeAIyJSiX9EzMfUvvtMEPac9v/et0Af38Yls7nUPUh1U7fTpZgIlxBMI1XdCewc9dzTIx5/A/hGaEsz8WBPTQtLctOZl5nidCkR66Yl/t7Lt2taKJmX4XA1JpLZnaLGMR6vj31nWtloZ+dXVTg7lYVZKXZh1EzIAt045khjJ90DHrsgOgERYePSOTYe3UzIAt04Zrj/fKMF+oQ2LplDa88gJy9aP7oZnwW6ccyemhZK5s4iLyPZ6VIi3vBfMTZ80VyNBbpxxKDHR/nZVhvdEqTC2WnkZ6fy1ulLTpdiIpgFunHEwbo2ege93LLsivvPzDhuXZbLW6dbbH50My4LdOOI3aeacbuEm+0MPWi3X5NHV7+HQw0dTpdiIpQFunHErpOXuKEom4wUm/88WLcsm4NLYNdJmwfJjM0C3cy4lu4BjjR2cHuJzeczGdlpSawqyGbXKQt0MzYLdDPj3qy+hCrcdo0F+mTdXpLLofp2OnqHnC7FRCALdDPjdp28RHZaIivzs5wuJercfk0ePoX/sdEuZgwW6GZGqSq7TzVzy7Jc3K6xpto3V7O6MJuM5AR2W7eLGYMFuplRJy520dQ1wB3Wfz4liW4XNy+bw66Tl2waAHMFC3Qzo4ZHaNx2jY0/n6rbSvI4195ny9KZK1igmxn1u+NNLJ+XwYIsW25uqu5c7v/r5rXjTQ5XYiKNBbqZMe29g+w/28bdpXOdLiWqFeSksWJBJq9WXXS6FBNhggp0EdkkIidEpFpEnhxj/8dF5HDg6y0RWR36Uk20e/1EM16fcveKeU6XEvXuXjGX8rOttPUMOl2KiSATBrqIuPEvK3cfUAo8IiKlo5qdAe5Q1VXAV4HtoS7URL9Xqy6SOyuZ1QXZTpcS9e5eMQ+fwmsnrNvFvCuYM/T1QLWq1qjqIPAcsGVkA1V9S1XbAptvAwWhLdNEu0GPjzdONHP3irm4bLjitK3Mz2JuRjK/sW4XM0IwgZ4P1I/Ybgg8N57PAL8aa4eIPC4i5SJS3txs42jjyb4zrXQPeKy7JURcLuGuFfPYdfISAx6v0+WYCBFMoI91OjXmAFgReR/+QP/SWPtVdbuqlqlqWV6ejUOOJ7+pukhygsumyw2he0rn0j3gYW9Nq9OlmAgRTKA3AIUjtguAxtGNRGQV8AywRVVtWRVzmary6rGL3FaSS2qS2+lyYsbNS3NJTXRbt4u5LJhA3w+UiMhiEUkCtgI7RjYQkSLgZ8AfqOrJ0JdpotnRxk7OtfdZd0uIpSS6uf2aXF45egGfz+4aNUEEuqp6gCeAV4Aq4HlVPSoi20RkW6DZV4A5wHdEpFJEysNWsYk6L71zHrdL+MB1850uJebcv3IBFzsHKK9tm7ixiXkJwTRS1Z3AzlHPPT3i8WPAY6EtzcQCVeWlw+e5ZVkuOelJTpcTc+5eMY+URBe/PNzI+sWznS7HOMzuFDVhdeRcJ3WtvTyw0s7OwyE9OYH3XzuXne9cwGvdLnHPAt2E1UvvnCfBJdxbaoEeLptXLeRS9wB7a2wsQryzQDdho6q89E6jdbeE2fuWzyUtyc0v3znvdCnGYRboJmwq6tupb+3jgVULnC4lpqUmublrxTxePnIBj9fndDnGQRboJmxeONBASqKL+6637pZw27xqAa09g7xZbUvTxTMLdBMW/UNeXjzUyAeum09GSqLT5cS8O5fnkZ2WyAsHzzldinGQBboJi98db6Kz38OHbrB52mZCcoKbLasX8srRC3T0DjldjnGIBboJixcONDA/M8XmbplBHykrZNDjY8fhK2bmMHHCAt2EXHPXAK+fbOahtfm4barcGXPdwkyunZ/BT8rrJ25sYpIFugm558vr8fqUj5RZd8tMEhE+WlbI4YYOTlzocroc4wALdBNSXp/yo7113LJsDkvzZjldTtx5aG0+SW4XP9xb63QpxgEW6CakXjvexLn2Pj6xodjpUuLS7PQkNq9ewAsHGujqt4uj8cYC3YTUD/bWMjcjmbtLbapcpzx60yJ6Br38zIYwxh0LdBMytS09vHGyma3ri0h027eWU1YXZrO6MJvv7zlr86THGfupMyHzzO4zJLpcfHxDkdOlxL1P3VxMTXMPu+3O0bhigW5CoqV7gOfL63l4bT7zMlOcLifu3b9yAXMzkvnuG6edLsXMoKACXUQ2icgJEakWkSfH2H+tiOwRkQER+ULoyzSR7vt7ahnw+Pij25c4XYrBf+foY7ct5q3TLVTU2WpG8WLCQBcRN/AUcB9QCjwiIqWjmrUCfwr8Q8grNBGvd9DDf+05yz2l81g214YqRorf31BMVmoi33ndztLjRTBn6OuBalWtUdVB4Dlgy8gGqtqkqvsBGycVh77/Vi3tvUNsu2Op06WYEWYlJ/DozYt49dhFTl60G43iQTCBng+MvJe4IfDcpInI4yJSLiLlzc3NUzmEiTCd/UM8/cZp7lyex7riHKfLMaN8+uZFpCe5+edXTzpdipkBwQT6WJNxTGkslKpuV9UyVS3Ly8ubyiFMhHlm9xk6+ob4wr3LnS7FjCEnPYk/un0JvzpywfrS40Awgd4AFI7YLgBsOjdDa88gz755hvuun8/1+VlOl2PG8dhtS5iTnsTXf3UcVRuXHsuCCfT9QImILBaRJGArsCO8ZZlo8A+/PkHfkJe/uPcap0sxVzErOYE/vauEvWdaee1Ek9PlmDCaMNBV1QM8AbwCVAHPq+pREdkmItsARGS+iDQAfw78tYg0iEhmOAs3znqnoYMf76vjUzcvYtncDKfLMRN4ZH0Ri3PT+dsXj9E/5HW6HBMmQY1DV9WdqnqNqi5V1f8TeO5pVX068PiCqhaoaqaqZgced4azcOMcn0/5yo4jzElP5vN3lzhdjglCUoKLr265ntqWXv7dhjHGLLtT1Ezaj/bVUVHXzpP3XUumrRcaNW4tyeXB1Qv599dPc+ZSj9PlmDCwQDeTUtvSw9d2VnFbSS4fumFKo1eNg/568wqSE1188SeH8Hh9TpdjQswC3QTN61O++JPDuEX4xodWIWLLy0WbuRkp/N1D11Ne22ZdLzHIAt0E7anXqtl3tpWv/F4pC7NTnS7HTNGWNfk8uHoh3/rtKRubHmMs0E1QXjvRxD//5iQPr83nw+tsrdBo99WHrmd+Zgp/8sODNHX1O12OCRELdDOhs5d6+PyPK7h2fiZfe3ildbXEgKzURLZ/ch1tvYP88Q8OMuCxoYyxwALdXFVTZz9/8Oxe3C7hu59YR2qS2+mSTIhctzCLf/jIag7UtvGXPz1sqxvFgASnCzCRq6N3iE8+u4+W7kF+9EcbKZqT5nRJJsQ2r1pIXWsv//flE6QlJfC1h6+3v8CimAW6GVNTVz+PPruf083dPPupG1lTmO10SSZM/uTOZfQMeHjqtdO4XfC3D16P22WhHo0s0M0Valt6ePTZfVzsHODZT93IbSU2M2as+8K9y/H4lO++UcOlrkG+tXUNKYnWvRZtrA/dvMdrx5v4vX97k7beIX7w2AYL8zghInz5vhX8zeZSXjl2gY9+dw/1rb1Ol2UmyQLdANA/5OVrO6v4w+/vpyAnjRefuNUWrIhDn7l1MU9/Yh1nLvVw/7/u5sVDjTblbhSxQDe8eeoS9//LbrbvqmHrjUW88Mc32wXQOPaB6+az809vY0luOp/7cQWPfb/cztajhDj127esrEzLy8sdeW/jd6i+nX969SRvnGymcHYqX//gKm5Zlut0WSZCeLw+vvfWWf7x1yfx+pSt6wv57PuWMS8zxenS4pqIHFDVsjH3WaDHl/4hL7+tauI//+cM5bVtZKYk8Ln3l/DJm4tJTrCLYOZKje19/NvvqvlJeT0ul7B55QJ+f0MR64pzbIijAyzQ41xH7xBvn2nh1WMXeeXIBboGPBTOTuXTNy/mozcWMivZBjuZidW19LJ992l+XtFI94CHxbnp3HvdPO4tncfqgmwS3NaDOxOmHegisgn4F8ANPKOqXx+1XwL77wd6gU+p6sGrHdMCPTz6h7ycvNhF1flOqs53cbCujSPnOvApZCQn8IHr5/Pg6oXcsizXxhqbKekZ8PDioUZeeuc8e0634PEpaUlu1hZls654NqULMlg2N4PiOWkkWsiH3LQCXUTcwEngHvwLRu8HHlHVYyPa3A98Dn+gbwD+RVU3XO24Fujj8/qUQY+PAY+XAY/v8uO+QR8dfUN09A3R3jdIR98QbT2DNHb0c66tj8b2Ppq7Bxj+X5qW5Ob6/CxuXjqHm5fmsqYwm6QE+wEzodPRN8TuU83sP9PK/rNtVF3ovPz9l+gWCnPSmJeZwvysFOZlpjA3I5nM1EQyUhLITPH/m5GSQFKCiyS3y/9v4LF154ztaoEezN/a64FqVa0JHOw5YAtwbESbLcB/qf+3w9siki0iC1T1/DRrv8IbJ5v56i+PXR5KdfnXkb7nnyv26+X9+t7tUb/PRr5u3Ncw+rXj7R+nhnFq96o/yD2TmFMjOcHFwuxUFmancOfyPBZmp7J8XgYrFmRSNDsNl52FmzDKSk1k86qFbF61EPCfvZ9u7qa6qZtTTd3UtvRwoaOffWdaudjZP6nv7eGAFwGXCK7AvzLisUv8Y+hl9HawbxJkw8n8FAXzi2jrjYU8dtuSSRw1OMEEej5QP2K7Af9Z+ERt8oH3BLqIPA48DlBUVDTZWgH/CubL5wUWJZb3/HP5P+S721ff/+7rZZz2Y+wb9eKxXnP1Y773f/bI9i6B5EQXSW534F8XyYkukhPcJCW4SElwkZWaSFZaItmpSWSnJdrdfCaipCcnsKogm1UF2Vfs8/mUjr4huvo9dPYP+b/6PHQPeBj0+Bj0eBn0+gKPfQx4fQx5FF/gBMinGvjynxT5fLxnW3l3OxjBXj+c1FXGIBvnzkqezFGDFkygj/XrZnTZwbRBVbcD28Hf5RLEe19hXXGO3fBiTBRyuYSc9CRy0pOcLiVmBdOh2gAUjtguABqn0MYYY0wYBRPo+4ESEVksIknAVmDHqDY7gE+K30agIxz958YYY8Y3YZeLqnpE5AngFfzDFp9V1aMisi2w/2lgJ/4RLtX4hy1+OnwlG2OMGUtQd5So6k78oT3yuadHPFbgs6EtzRhjzGTYoGRjjIkRFujGGBMjLNCNMSZGWKAbY0yMcGy2RRFpBmqn+PJc4FIIy3GSfZbIZJ8l8sTK54DpfZZiVR1zbUjHAn06RKR8vMlpoo19lshknyXyxMrngPB9FutyMcaYGGGBbowxMSJaA3270wWEkH2WyGSfJfLEyueAMH2WqOxDN8YYc6VoPUM3xhgzigW6McbEiKgNdBH5pogcF5HDIvLfIpLtdE1TJSIfEZGjIuITkagbliUim0TkhIhUi8iTTtczHSLyrIg0icgRp2uZDhEpFJHXRKQq8L31eadrmioRSRGRfSJyKPBZ/tbpmqZLRNwiUiEivwzlcaM20IFXgetVdRX+Ray/7HA903EE+CCwy+lCJiuwiPhTwH1AKfCIiJQ6W9W0fA/Y5HQRIeAB/kJVVwAbgc9G8f+XAeD9qroaWANsCqy7EM0+D1SF+qBRG+iq+mtV9QQ238a/SlJUUtUqVT3hdB1TdHkRcVUdBIYXEY9KqroLaHW6julS1fOqejDwuAt/eOQ7W9XUqF93YDMx8BW1ozlEpAB4AHgm1MeO2kAf5Q+BXzldRJwab4FwEyFEZBGwFtjrcClTFuiiqASagFdVNWo/C/At4C8BX6gPHNQCF04Rkd8A88fY9Veq+otAm7/C/+flD2eytskK5rNEqaAWCDfOEJFZwAvAn6lqp9P1TJWqeoE1gWtl/y0i16tq1F3nEJHNQJOqHhCRO0N9/IgOdFW9+2r7ReRRYDNwl0b4gPqJPksUswXCI5SIJOIP8x+q6s+cricUVLVdRF7Hf50j6gIduAV4UETuB1KATBH5gap+IhQHj9ouFxHZBHwJeFBVe52uJ44Fs4i4mWEiIsB/AFWq+k9O1zMdIpI3PIpNRFKBu4HjjhY1Rar6ZVUtUNVF+H9WfheqMIcoDnTg20AG8KqIVIrI0xO9IFKJyMMi0gDcBLwkIq84XVOwAhemhxcRrwKeV9WjzlY1dSLyY2APsFxEGkTkM07XNEW3AH8AvD/w81EZOCuMRguA10TkMP4TiFdVNaTD/WKF3fpvjDExIprP0I0xxoxggW6MMTHCAt0YY2KEBboxxsQIC3RjjIkRFujGGBMjLNCNMSZG/H9050A72M9pqgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "手工计算分布函数的值： 0.8413447460685428\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate as integ\n",
    "\n",
    "# 正态分布非常重要，很多统计分析和机器学习的算法都使用该分布\n",
    "\n",
    "# 构建服从正态分布的概率密度函数的值，均值mu为1，标准差sigma为0.5\n",
    "x = np.arange(-2, 4, 0.01)\n",
    "# 根据x计算正态分布的概率密度函数\n",
    "n_density = (1/(((2*np.pi)**0.5)*0.5))*np.exp((-(x-1)**2)/(2*(0.5**2)))  # 第一个公式f(x)\n",
    "# 正态分布的密度函数曲线，后面还有更具体的图形例子\n",
    "print('正态分布概率密度函数曲线（均值：1，标准差：0.5）：')\n",
    "plt.plot(x, n_density)\n",
    "plt.show()\n",
    "\n",
    "# 累积分布函数的计算，设x0 = 1.5，则F(x0) = P(x<1.5)\n",
    "# 即求概率密度函数f(x) 在 (-Inf, 1.5) 区间的积分\n",
    "# 被积分函数\n",
    "def nf(x, mu, sigma):\n",
    "    return (1/(((2*np.pi)**(0.5))*sigma))*np.exp((-(x-mu)**2)/(2*(sigma**2)))\n",
    "\n",
    "# 累积分布函数是对概率密度函数求积分（针对连续型随机变量，离散型随机变量是求和）\n",
    "# nf 概率密度函数\n",
    "# float('-Inf'), 1.5  -> P(x<1.5) ，是nf函数的参数x\n",
    "# args=(1, 0.5)  -> 1 是mu，0.5是sigma\n",
    "n_distrib = integ.quad(nf, float('-Inf'), 1.5, args=(1, 0.5))[0]  # 第二个公式F(x)的实现\n",
    "\n",
    "print('\\n手工计算分布函数的值：', n_distrib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6fe4693e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Scipy函数计算正态分布的累积分布函数的值： 0.8413447460685429\n",
      "\n",
      "服从正态分布的随机数生成： [0.89303624 0.47842399 0.78380919 1.21025428 0.63242221]\n",
      "\n",
      "概率密度函数（显示前10个）： [1.21517657e-08 1.36983376e-08 1.54355684e-08 1.73861599e-08\n",
      " 1.95754158e-08 2.20315272e-08 2.47858886e-08 2.78734461e-08\n",
      " 3.13330824e-08 3.52080407e-08]\n",
      "\n",
      "生存函数： 0.15865525393145707\n",
      "\n",
      "分位数函数： 1.4972289416048765\n",
      "\n",
      "置信度0.95的置信区间： (0.020018007729972975, 1.979981992270027)\n"
     ]
    }
   ],
   "source": [
    "# 通过scipy计算\n",
    "\n",
    "# scipy的子包stats下的norm表示正态分布，norm包含几乎所有正态分布计算的函数\n",
    "# 类似的还有t、f、chi2等子包对应其他各种分布\n",
    "\n",
    "# 直接调用Scipy的累积分布函数\n",
    "# 1  -> 均值\n",
    "# 0.5  -> 标准差\n",
    "n_distrib_scipy = st.norm.cdf(1.5, 1, 0.5)\n",
    "print('\\nScipy函数计算正态分布的累积分布函数的值：', n_distrib_scipy)\n",
    "\n",
    "# Scipy有关概率分布的其他应用\n",
    "# 1 均值，0.5 标准差，5 生成5个\n",
    "print('\\n服从正态分布的随机数生成：', st.norm.rvs(1, 0.5, 5))\n",
    "print('\\n概率密度函数（显示前10个）：', st.norm.pdf(x, 1, 0.5)[0: 10])\n",
    "print('\\n生存函数：', st.norm.sf(1.5, 1, 0.5))  # sf = 1 - cdf\n",
    "print('\\n分位数函数：', st.norm.ppf(0.84, 1, 0.5))  # 注意：通常都是通过下分位数（下尾）计算分位数。\n",
    "print('\\n置信度0.95的置信区间：', st.norm.interval(0.95, 1, 0.5))  # 该分布下随机变量有95%的可能性位于该区间\n",
    "# 均值为1，标准差为0.5的正态分布生成的随机变量，有95%的概率落在置信区间"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4057fd3a",
   "metadata": {},
   "source": [
    "#### T分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e8d59149",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFlCAYAAAAki6s3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7I0lEQVR4nO3deXyV9Z3+/9c7JxuEhCUrhAABEiCsIoKCgKIirlRrW62dtmOt41i/nU6nv9ZOO52l00477bSdtlprW506rbWO1pYqCirK4oIERNaQhLAkQJKTQEhCyHo+vz8SaMQgJ5DkPsv1fDx4mHPOfedct5Bz5f7c9/25zTmHiIiIhKYYrwOIiIjI2amoRUREQpiKWkREJISpqEVEREKYilpERCSEqahFRERCWKzXAXqSlpbmxo0b53UMERGRAbF58+Ya51x6T6+FZFGPGzeOwsJCr2OIiIgMCDM7cLbXNPQtIiISwlTUIiIiIUxFLSIiEsJU1CIiIiFMRS0iIhLCVNQiIiIhTEUtIiISwlTUIiIiIUxFLSIiEsJU1CIiIiFMRS0iIhLCQnKubxEZWDWNLeytbsQB49OSyEhJ9DqSiHQJqqjNbBnw34AP+KVz7jtnvL4c+CYQANqBLzjnNnS9th9oADqAdufcnD5LLyIXZG2xn5+uKWHT/mPveX5WzjDuu2IC1xRkYmYepRMRCKKozcwHPAhcA1QAm8xshXNuV7fFXgFWOOecmc0AngImd3v9SudcTR/mFpEL0NzWwT8+u50/bDnE6OGD+NLSfGbmDCPGjG0Vx3mqsJx7/ncz10/P4nu3zSQpQYNvIl4J5qdvLlDqnCsDMLMngeXA6aJ2zjV2Wz4JcH0ZUkT6Tl1TK3f9zya2HKzj81flcf+VE4mP/cvpKgsmpvHZhbk8sr6M76/aw4HaJh7760vISNZwuIgXgjmZLBso7/a4ouu59zCzW8ysCHgeuKvbSw5YbWabzeyeCwkrIhemua2Dz/y6kB2H6nnoztl88Zr895T0KbG+GO67YiK/+tQllPlP8OlHN9HQ3OZBYhEJpqh7OkD1vj1m59yzzrnJwIfoPF59ygLn3GzgOuBzZraoxzcxu8fMCs2s0O/3BxFLRHrDOcff/34rWw4e40e3z+L66SPPuc6VkzN46BOz2VPVwH2/3UJHQINlIgMtmKKuAHK6PR4NHD7bws65dcAEM0vreny467/VwLN0DqX3tN4jzrk5zrk56enpQcYXkWD9+o39vLCjkn+8bkpQJX3KlZMy+ObyaawvqeFnr5X2Y0IR6UkwRb0JyDOzXDOLB24HVnRfwMwmWtepoWY2G4gHas0sycySu55PApYCO/pyA0Tk3Ioq6/n2yiKumpzB3Qtze73+HXNzuHnmKH74cglbDh479woi0mfOWdTOuXbgfmAVsBt4yjm308zuNbN7uxb7MLDDzLbSeYb4x5xzDsgENpjZu8DbwPPOuRf7YTtE5Cw6Ao4HntlOcmIs3/vIzPO63MrM+NYt08hKSeSrz2ynrSPQD0lFpCdBXXPhnFsJrDzjuYe7ff1d4Ls9rFcGzLzAjCJyAZ54+yBby+v44cdmMiIp/ry/T3JiHP9y81Q++3ghv9qwj3sXT+jDlCJyNppCVCSCHW9q43svFjF/QiofmvW+izV67ZqCTJYWZPLfL5dQXd/cBwlF5FxU1CIR7Ofr9lLf3M7XbyjosxnGvnbDFNo6AvzolZI++X4i8sFU1CIRqrq+mUdf38fyWaMoGJXSZ993bGoSd84bw+83lbPX33juFUTkgqioRSLUj9eU0N7h+OI1+X3+vf/fVXkkxMbwE+1Vi/Q7FbVIBKo83szvN5XzsUtyGJua1OffP21IAnfOG8Oftx3hYG1Tn39/EfkLFbVIBHrsjX10BFy/npl998Lx+Mz4+bq9/fYeIqKiFok49c1tPPHWQW6YMYqcEYP77X0yUxK5bc5o/q+wQmeAi/QjFbVIhPndxoM0tLTzN4vG9/t73btoAu2BAI+9sb/f30skWqmoRSJIW0eAR1/fx4KJqUzLHtrv7zcmdTDXFGTy5NsHaW7r6Pf3E4lGKmqRCPLyriqq6lv46/m9n8/7fH3qsnEca2rj+W1HBuw9RaKJilokgvxm4wGyhw3iyskZA/ael01IZWLGEB5/c/+AvadINFFRi0SIMn8jr5fWcsfcHHwxfTMLWTDMjE9eNpZ3K46ztbxuwN5XJFqoqEUixBMbDxIbY3z0kpxzL9zHbp09mqR4H79968CAv7dIpFNRi0SAlvYOnt5SwbVTs8hIThzw9x+SEMsNM0by/PYjnGhpH/D3F4lkKmqRCLBmdzV1TW18zIO96VNuuziHptYOXtxR6VkGkUikohaJAM9sqSAjOYEFE9M8y3DJuOGMTR3M05srPMsgEolU1CJhrqaxhdf2+LnlouwBPYnsTGbGbbNH82ZZLeVHNf+3SF9RUYuEuRVbD9MecNw6e7TXUbj14tGYde7hi0jfUFGLhLk/vFPBtOwUJmUlex2F7GGDWDAhjac3VxAIOK/jiEQEFbVIGNtT2cCOQ/XcepH3e9OnfPjibCqOneSd8mNeRxGJCCpqkTD27DuH8MUYN88a5XWU066ekklCbAx/fldTior0BRW1SJhyzvH89sMsmJhG2pAEr+OclpwYx5WTMnh++xE6NPwtcsFU1CJhasehesqPnuSG6VleR3mfm2aOwt/QwsZ9tV5HEQl7KmqRMPX89iPExhhLC0KvqJdMzmBwvI/ndEctkQumohYJQ6eGvedPTGN4UrzXcd5nULyPq6dk8sL2I7R1BLyOIxLWVNQiYejUsPeN00d6HeWsbpwxkmNNbbyxV8PfIhdCRS0Shk4Pe0/N9DrKWS2elE5yQizPvXvY6ygiYU1FLRJmup/tPWxw6A17n5IQ62PJlAxe3l1Fu4a/Rc6bilokzPzlbO/QHfY+5dqpWRxraqPwgCY/ETlfKmqRMLNyR+gPe5+yOD+d+NgYVu3UrS9FzpeKWiTMrN5ZyaXjU0N62PuUpIRYFk5MY/XOKpzT5Cci50NFLRJGyvyN7PWf4JqC0N+bPuXaqVkcqjvJzsP1XkcRCUsqapEw8tKuKgCumpLhcZLgXTUlgxjrHAkQkd5TUYuEkZd3V1EwMoXRwwd7HSVoqUMSmDNuBKu7fskQkd5RUYuEidrGFjYfOBZWw96nXDs1i6LKBg7UnvA6ikjYCaqozWyZme0xs1Ize6CH15eb2TYz22pmhWZ2ebDrikhwXimqJuAIy6Je2pVZZ3+L9N45i9rMfMCDwHVAAXCHmRWcsdgrwEzn3CzgLuCXvVhXRILw0q4qRg1NZOqoFK+j9FrOiMEUjEzh5V3VXkcRCTvB7FHPBUqdc2XOuVbgSWB59wWcc43uL9deJAEu2HVF5NxOtnawvsTP1QWZmJnXcc7LkskZbD54jONNbV5HEQkrwRR1NlDe7XFF13PvYWa3mFkR8Dyde9VBr9u1/j1dw+aFfr8/mOwiUeP10hqa2wJhOex9ypIpGXQEHGtL9PMt0hvBFHVPv76/b+YC59yzzrnJwIeAb/Zm3a71H3HOzXHOzUlPTw8ilkj0eHl3FUMSYpmXm+p1lPM2c/QwRiTFs2a3zv4W6Y1giroCyOn2eDRw1tvhOOfWARPMLK2364rI+znnWFvs5/KJacTHhu+FGr4Y44r8dNYW++kIaJYykWAF81O/Ccgzs1wziwduB1Z0X8DMJlrXgTMzmw3EA7XBrCsiH6y0upEjx5tZlB/+I01LpmRwrKmNreW6SYdIsGLPtYBzrt3M7gdWAT7gUefcTjO7t+v1h4EPA580szbgJPCxrpPLely3n7ZFJCKtLe48prsoP83jJBduYV46vhjjld3VXDx2hNdxRMLCOYsawDm3Elh5xnMPd/v6u8B3g11XRIK3rqSG8elJYTUb2dkMHRTHnLHDWVNUzZeXTfY6jkhYCN8DXiJRoLmtg41ltSzKC/9h71OumpJBUWUDh+tOeh1FJCyoqEVC2Nv7jtLSHmBxBByfPmXJ5M4biqwp0uQnIsFQUYuEsHXFfuJ9McwbHznHcyekDyFnxCBeVVGLBEVFLRLC1pX4uSR3OIPjgzqdJCyYGUsmZfD63hqa2zq8jiMS8lTUIiHqyPGTFFc1RtTx6VOunJxBc1uAN8tqvY4iEvJU1CIhan1xDUBEXD99pkvHp5IQG8O6Yk0nKnIuKmqRELW2xE9GcgKTs5K9jtLnEuN8zBufqqIWCYKKWiQEdQQcG0pqWJiXHrZ3yzqXxfnp7PWfoOJYk9dRREKailokBG2rqOP4ybaImI3sbBZ3bdu6riF+EemZilokBK0rrsGsc8rNSDUhfQjZwwaxtliXaYl8EBW1SAhaV+JnevZQRiTFex2l35gZi/LTeKO0lraOgNdxREKWilokxBw/2cbW8rqImo3sbBbnp9PQ0s47B+u8jiISslTUIiHmjdIaOgIuIi/LOtP8iWn4Ykxnf4t8ABW1SIhZV+InOSGWWTnDvI7S71IS45g9ZtjpW3mKyPupqEVCiHOOdcU1zJ+YSpwvOn48F+Wls/3QcWoaW7yOIhKSouOTQCRM7PWf4FDdyagY9j5l8aTObd1Qosu0RHqiohYJIaeO1Ubi/N5nM21U59ntGv4W6ZmKWiSErCvxMz4tiZwRg72OMmBiYoyFeWmsL/ETCDiv44iEHBW1SIhobuvgrbLaqBr2PmVxfjo1ja3sOlLvdRSRkKOiFgkRhfuP0dwWiOhpQ8/m1AxsGv4WeT8VtUiIWFfiJ94Xw6XjU72OMuDSkxOYOipFRS3SAxW1SIhYV+xnzrjhDI6P9TqKJxblp7PlwDEamtu8jiISUlTUIiGgqr6ZosqGqDw+fcrCvDTaA463yo56HUUkpKioRUJANF6WdaY5Y0cwON6n6URFzqCiFgkB60pqSE9OYMrIZK+jeCY+NobLxqeyrkRFLdKdilrEYx0Bx/oSPwvz0jAzr+N4amFeGgdqmzhQe8LrKCIhQ0Ut4rHth45T19QWFbe1PJdTx+jXaTpRkdNU1CIeW1fsxwwunxh910+fKTctidHDB+k4tUg3KmoRj60r9jNt1FBShyR4HcVzZsbCvHTe3FtLW0fA6zgiIUFFLeKh+uY23imvi8rZyM5mcX4ajS3tvHOwzusoIiFBRS3ioTdKa+gIOBbnZ3gdJWRcNiENX4xp+Fuki4paxENri2sYkhDLRWOGeR0lZAwdFMesnGG6TEuki4paxCPOOdYV+5k/IZU4n34Uu1uUl872Q8c5eqLV6ygintOng4hHympOcKjuZFRPG3o2i/LTcA42lOoyLZGgitrMlpnZHjMrNbMHenj9TjPb1vXnDTOb2e21/Wa23cy2mllhX4YXCWenjsHq+un3mzF6GEMHxek4tQhwztv0mJkPeBC4BqgANpnZCufcrm6L7QMWO+eOmdl1wCPAvG6vX+mc06/GIt2sK/aTm5ZEzojBXkcJOb4Y4/KJaawv8eOci/oZ2yS6BbNHPRcodc6VOedagSeB5d0XcM694Zw71vXwLWB038YUiSwt7R28VXaURXm6LOtsFuWnUVXfQnFVo9dRRDwVTFFnA+XdHld0PXc2nwFe6PbYAavNbLOZ3dP7iCKRp3D/MU62dej49AdY2HUnMQ1/S7QLpqh7GnNyPS5odiWdRf2Vbk8vcM7NBq4DPmdmi86y7j1mVmhmhX6/fjAlsq0r9hPnMy4dn+p1lJA1atggJmYM0WVaEvWCKeoKIKfb49HA4TMXMrMZwC+B5c652lPPO+cOd/23GniWzqH093HOPeKcm+Ocm5Oerr0MiWxri/3MGTuCpIRzniYS1RblpfP2vqM0t3V4HUXEM8EU9SYgz8xyzSweuB1Y0X0BMxsD/AH4K+dccbfnk8ws+dTXwFJgR1+FFwlH1fXNFFU2aNg7CIvy02hpD7Bx31Gvo4h45py/zjvn2s3sfmAV4AMedc7tNLN7u15/GPgGkAo81HV2Zrtzbg6QCTzb9Vws8IRz7sV+2RKRMHHqFo6a3/vc5uWmEh8bw/pivy5jk6gV1Libc24lsPKM5x7u9vXdwN09rFcGzDzzeZFotq7YT9qQBKZkpXgdJeQNivcxd9wIHaeWqKaZyUQGUCDg2FBaw6K8NGJidG1wMBblp1Fc1ciR4ye9jiLiCRW1yADacbhz/modnw7eqcu01pdoziSJTipqkQF06prgyzXRSdAmZyWTkZyg66klaqmoRQbQuuIapmWnkDYkwesoYcPMWJiXzoaue3eLRBsVtcgAaWhuY8vBYyzK07B3by3KT6OuqY3th457HUVkwKmoRQbIG3traQ84HZ8+D5dPTMMM1mv4W6KQilpkgKwr9jMkIZbZY4Z7HSXspA5JYNqoobpMS6KSilpkADjnWFfi57IJnRN4SO8tyk9jy8E66pvbvI4iMqD0iSEyAPbXNlF+9KSGvS/Awrx0OgKON/fWnnthkQiiohYZAKcuLVqsE8nO2+wxw0mK9+kyLYk6KmqRAbCu2M+41MGMSR3sdZSwFR8bw2UT0lhX4sc5XaYl0UNFLdLPWto7eLOsVsPefWBRfhrlR09yoLbJ6ygiA0ZFLdLPNu8/RlNrh+7+1AdOXYOus78lmqioRfrZ2mI/8b4YLh2f6nWUsDcuLYkxIwbrOLVEFRW1SD9bW+znktzhJCUEdVdZOYeFeWm8ubeW1vaA11FEBoSKWqQfVR5vpqiyQcPefWhRfjonWjvYcvCY11FEBoSKWqQfnb4sKz/D4ySRY/6EVGJjTMPfEjVU1CL9aG2xn6yURPIzh3gdJWIkJ8Yxe8xwnVAmUUNFLdJP2jsCrC/xszg/HTPzOk5EWZiXxo5D9dQ2tngdRaTfqahF+sm7FXXUN7ezeJKOT/e1U9ekbyit8TiJSP9TUYv0k7V7/PhijAUT07yOEnGmZQ9l+OA41uo4tUQBFbVIP1lb7GdWzjCGDorzOkrEOfUL0PqSGk0nKhFPRS3SD2obW9h26Lguy+pHi/LT8Te0UFTZ4HUUkX6lohbpBxtKa3AOFXU/Oj2dqIa/JcKpqEX6wdo9fkYkxTM9e6jXUSJW1tDOy97Wl+iEMolsKmqRPhYIONaV+FmYl0ZMjC7L6k+L8tJ5e/9RTrZ2eB1FpN+oqEX62K4j9dQ0tmrYewAsyk+ntT3AW/tqvY4i0m9U1CJ97NQlQwvzVNT9bW7uCBJiY1hfrOFviVwqapE+tnaPn2nZKaQnJ3gdJeIlxvmYmztC04lKRFNRi/Sh4yfb2HzwmIa9B9Di/HRKqxs5XHfS6ygi/UJFLdKH1hX76Qg4lkzW3bIGyqnpRHWZlkQqFbVIH3q1qJrhg+OYlTPc6yhRIy9jCCOHJvLaHhW1RCYVtUgf6Qg4Xt1TzRWTMvDpsqwBY2YsmZzB+hI/Le26TEsij4papI9sLa/jWFMbV2rYe8BdNSWDE60dvL3vqNdRRPqcilqkj6wpqsIXYyzWZVkDbv6ENBLjYnhld7XXUUT6XFBFbWbLzGyPmZWa2QM9vH6nmW3r+vOGmc0Mdl2RSLGmyM/FY4czdLDuljXQEuN8zJ+QxitFVbqblkSccxa1mfmAB4HrgALgDjMrOGOxfcBi59wM4JvAI71YVyTsHa47ye4j9VylYW/PLJmcQfnRk+z1N3odRaRPBbNHPRcodc6VOedagSeB5d0XcM694Zw71vXwLWB0sOuKRIJX93QOueqyLO+c+n+v4W+JNMEUdTZQ3u1xRddzZ/MZ4IXermtm95hZoZkV+v26zELCy5rd1eSMGMTEjCFeR4lao4YNYsrIFF4pUlFLZAmmqHu6zqTHg0BmdiWdRf2V3q7rnHvEOTfHOTcnPV0n40j4aG7r4PW9NSyZlIGZLsvy0lWTM9h84Bh1Ta1eRxHpM8EUdQWQ0+3xaODwmQuZ2Qzgl8By51xtb9YVCWdv7q2luS3AkimZXkeJekumZNARcKdvjCISCYIp6k1Anpnlmlk8cDuwovsCZjYG+APwV8654t6sKxLu1hRVMyjOx7zcEV5HiXozRw8jNSmeNRr+lggSe64FnHPtZnY/sArwAY8653aa2b1drz8MfANIBR7qGvpr7xrG7nHdftoWkQHnnGNNUTULJqaRGOfzOk7U88UYV0zK4OXdVbR3BIj1aaoICX/nLGoA59xKYOUZzz3c7eu7gbuDXVckUhRXNXKo7iT3L5nodRTpctWUDJ7ZUsGWg3XM1SiHRAD9uilyAV4pqgLgykm6LCtULMxLIzbGTv/diIQ7FbXIBXhpVxXTs4eSNTTR6yjSJTkxjrm5I1ij66klQqioRc5TdX0z7xysY2mBzvYONUsmZ1BS3ciB2hNeRxG5YCpqkfP00u7OodWlU7M8TiJnWlrQ+Xfy0i4Nf0v4U1GLnKeXdlUxZsRg8jM1G1moGZM6mMlZyazaWel1FJELpqIWOQ8NzW28UVrL0oJMzUYWopZOzaLwwDFqGlu8jiJyQVTUIudhbbGf1o6Ahr1D2LVTM3EOXtmt4W8JbypqkfOwemcVI5LiuXjscK+jyFkUjEwhe9ggVu1UUUt4U1GL9FJre4BX91Rz9ZQMfDEa9g5VZsa1U7PYUFpDY0u713FEzpuKWqSXNu6rpaG5nWsKNOwd6pZOzaS1PcA63aRDwpiKWqSXVu+sYlCcj4V5aV5HkXOYM3Y4wwfH6exvCWsqapFeCAQcL+2qYlG+bsIRDmJ9MVw9JZM1RdW0tge8jiNyXlTUIr2w/dBxKuubT0+oIaFv6dQsGprbeaus1usoIudFRS3SC6t3VeKLMZZM1k04wsXCvDQGxflYvUvD3xKeVNQiQXLO8cL2SubljmB4UrzXcSRIiXE+Fuen89KuKgIB53UckV5TUYsEaU9VA2U1J7h++kivo0gvLZ2aSVV9C++U13kdRaTXVNQiQVq57QgxBtdqNrKwc9WUTOJ9MazcfsTrKCK9pqIWCdLKHZXMzR1BenKC11Gkl4YOimNRfhortx/R8LeEHRW1SBCKqxoorW7kBg17h63rp4/kyPFmDX9L2FFRiwRh5fYjmMG10zTsHa6uLugc/n5+m4a/JbyoqEWCsHL7ES4ZN4KM5ESvo8h5SkmMY1F+Oi/s0PC3hBcVtcg5lFY3UFylYe9IcMOMrK7h72NeRxEJmopa5BxWbq/EDJZp2DvsXT0lk/jYGJ7T8LeEERW1yDms3H6EOWOHk5miYe9wl5wYx+L8dF7YXqnhbwkbKmqRD7DX30hRZQPXTdOwd6S4ccZIKuub2XJQw98SHlTUIh9gxdbDmHV+uEtkuErD3xJmVNQiZ+GcY8W7h5k/IZUMDXtHjCEJsVyhs78ljKioRc5i+6Hj7Ks5wc0zR3kdRfrYDTNGUlXfwqb9R72OInJOKmqRs/jT1sPE+2JYNlXD3pHm6imZDIrz8ceth72OInJOKmqRHnQEHH9+9zBXTEpn6OA4r+NIH0tKiOXaqZms3H6E1vaA13FEPpCKWqQHG8tqqW5oYfmsbK+jSD9ZflE2x0+28dqeaq+jiHwgFbVID1a8e5ikeB9XTcnwOor0k4UT00hNiuePWw95HUXkA6moRc7Q0t7Byu1HuHZqFolxPq/jSD+J9cVw08xRvLy7mvrmNq/jiJyVilrkDGv3+KlvbufmWTrbO9ItnzWK1vYAL26v9DqKyFkFVdRmtszM9phZqZk90MPrk83sTTNrMbMvnfHafjPbbmZbzaywr4KL9Jc/bT3MiKR4FkxM8zqK9LNZOcMYlzpYw98S0s5Z1GbmAx4ErgMKgDvMrOCMxY4Cnwe+f5Zvc6VzbpZzbs6FhBXpb8eb2nhpVxU3zxxFnE8DTpHOzFg+K5s3y2qpPN7sdRyRHgXzSTQXKHXOlTnnWoEngeXdF3DOVTvnNgE60CNhbcW2w7R2BLjt4tFeR5EB8qGLsnEOVryrvWoJTcEUdTZQ3u1xRddzwXLAajPbbGb39CacyEB7enMFk7OSmToqxesoMkBy05KYmTOMP2w5hHOaUlRCTzBFbT0815t/zQucc7PpHDr/nJkt6vFNzO4xs0IzK/T7/b349iJ9o7S6gXfL67jt4tGY9fTPXiLVbRePpqiygR2H6r2OIvI+wRR1BZDT7fFoIOh595xzh7v+Ww08S+dQek/LPeKcm+Ocm5Oenh7stxfpM09vPoQvxjTJSRS6eeYoEmJjeKqw/NwLiwywYIp6E5BnZrlmFg/cDqwI5pubWZKZJZ/6GlgK7DjfsCL9pSPgePadCq6clE56coLXcWSADR0Ux7VTs/jT1kM0t3V4HUfkPc5Z1M65duB+YBWwG3jKObfTzO41s3sBzCzLzCqALwJfN7MKM0sBMoENZvYu8DbwvHPuxf7aGJHztb7ET1V9i04ii2IfnZNDfXM7q3dVeR1F5D1ig1nIObcSWHnGcw93+7qSziHxM9UDMy8koMhAeHpzBcMHx7FkcqbXUcQj8yekkj1sEP9XWK5bm0pI0YWiEvWON7WxelcVy2dlEx+rH4loFRNj3HbxaDaU1nCo7qTXcURO06eSRL0/bj1Ea7uunZbOs7+dg2c2V3gdReQ0FbVENeccv3v7INOzhzIte6jXccRjOSMGM39CKk9vriAQ0DXVEhpU1BLVthyso6iygY/PG+N1FAkRH52Tw8GjTby1r9brKCKAilqi3O/ePkhSvE8nD8lpy6ZlMXRQHL/deNDrKCKAilqi2PGTbTy37TDLL8omKSGoCyAkCiTG+bjt4tGs2lFJdYNu1CHeU1FL1Hp2SwXNbQE+PlfD3vJed84bQ3vA8dQmzVQm3lNRS1TqPImsnJmjdRKZvN/49CEsmJjK794up0MnlYnHVNQSlbYcPMaeqgbu0N60nMWd88ZyqO4kr+2p9jqKRDkVtUSl3248yJCEWG7SSWRyFtcUZJKRnMBv3jrgdRSJcipqiTq1jS08t+0It+gkMvkAcb4Ybr8kh9eK/ZQfbfI6jkQxFbVEnd+9fZDW9gCfmj/W6ygS4m6fOwYDnnhbl2qJd1TUElXaOgL871sHWJiXxsSMZK/jSIgbNWwQSyZn8vtN5br9pXhGRS1R5cUdlVTVt3DXglyvo0iYuGvBOI6eaGXF1sNeR5EopaKWqPLY6/sYlzqYxfnpXkeRMHHZhFQmZyXz6Ov7cE6XasnAU1FL1Hi3vI4tB+v41PxxxMSY13EkTJgZd12eS1FlA6+Xav5vGXgqaokav35jP0MSYnU7S+m1m2eOIm1IPI++vs/rKBKFVNQSFaobmvnztsPcdvFokhPjvI4jYSYxzsed88aypqiavf5Gr+NIlFFRS1R4/I0DtAccn5o/zusoEqY+celY4n0xPKa9ahlgKmqJeCda2nn8zf1cW5BFblqS13EkTKUnJ3DzrFE8s/kQdU2tXseRKKKiloj3u7cPUt/czt8sHu91FAlzn7k8l5NtHbpXtQwoFbVEtNb2AL/asI95uSO4aMxwr+NImJsyMoVF+ek8umEfJ1s1AYoMDBW1RLQV7x7myPFm7r1igtdRJEJ87ooJ1J5o5alC3ataBoaKWiJWIOD4+dq9TM5K5gpNcCJ9ZG7uCOaMHc4j68po6wh4HUeigIpaItaaompKqhu5d/EEzDTBifQNM+O+KydwqO4kf9K0ojIAVNQSkZxzPLx2L9nDBnHDjJFex5EIc+WkDCZnJfOz10oJBDStqPQvFbVEpDf21lJ44Bj3LBpPnE//zKVvde5VT2Sv/wSrd1V6HUcinD7BJOI45/jRy8VkpSTysUtyvI4jEeqG6SMZlzqYB1/dq5t1SL9SUUvEeb20lk37j3HflRNIjPN5HUcilC/G+NsrJrD90HHWFFV7HUcimIpaIor2pmUg3Tp7NGNTB/Nfq4t1rFr6jYpaIsqG0hoKDxzjc1dOICFWe9PSv+J8MfzdVXnsOlLPqp06Vi39Q0UtEaNzb7qEkUMT+aj2pmWALJ+VzYT0JH74cjEd2quWfqCiloixvqSGzQeOcd+VE7U3LQPGF2N84ep8iqsaeW6brquWvqeilogQCDj+c1UR2cMG8dE5o72OI1HmhukjmZyVzI9eLqFds5VJH1NRS0T487bD7DhUzz8szdfetAy4mBjj76/JZ1/NCf7wziGv40iECaqozWyZme0xs1Ize6CH1yeb2Ztm1mJmX+rNuiIXqqW9g++v3sOUkSl8aFa213EkSi0tyGTm6KH8YHWx7qwlfeqcRW1mPuBB4DqgALjDzArOWOwo8Hng++exrsgF+e1bByk/epKvLJtETIzm9BZvmBn/eP0UKuub+dWGMq/jSAQJZo96LlDqnCtzzrUCTwLLuy/gnKt2zm0C2nq7rsiFqG9u4ydrSpg/IZXFukOWeGze+FSWFmTys9f24m9o8TqORIhgijob6H7j1Yqu54IR9Lpmdo+ZFZpZod/vD/LbS7R7ZG0Zx5ra+Op1U3SHLAkJD1w3mZb2AD96udjrKBIhginqnj79gr1YMOh1nXOPOOfmOOfmpKdrz0jO7XDdSX65oYybZo5i+uihXscRAWB8+hDunDeGJzeVU1LV4HUciQDBFHUF0H32iNFAsBcLXsi6Ih/o2yt34xx8+dpJXkcReY/PX5XH4Dgf33mhyOsoEgGCKepNQJ6Z5ZpZPHA7sCLI738h64qc1Zt7a3lu2xHuXTyBnBGDvY4j8h6pQxK478qJvFJUzdpiHcqTC3POonbOtQP3A6uA3cBTzrmdZnavmd0LYGZZZlYBfBH4uplVmFnK2dbtr42R6NDeEeBf/7yT7GGD+NsrJngdR6RHd10+jty0JP5lxU5a2nW5lpy/2GAWcs6tBFae8dzD3b6upHNYO6h1RS7EbzcepKiygZ/dOVu3sZSQlRDr419unsqnHn2bX67fx+eunOh1JAlTmplMwsrRE6381+o9LJiYyrJpWV7HEflAi/PTuW5aFj9ZU0LFsSav40iYUlFLWPnuC0WcaO3gn2+aqsuxJCx8/cYCDOObz+3yOoqEKRW1hI0399by+8Jy7r48l/zMZK/jiAQle9gg7l8ykVU7q3htT7XXcSQMqaglLDS3dfC1Z7czZsRgvnB1vtdxRHrl7oW5jE9P4p/+tIOm1nav40iYUVFLWPjpmlLKak7wrVumMSheJ5BJeEmI9fGdW2dQfvQk31u1x+s4EmZU1BLyiirreXjtXm6dnc3CPM1aJ+Fpbu4IPnnZWP7njf1sPnDU6zgSRlTUEtI6Ao4HntlOyqA4vn6Dbrwm4e3LyyYzauggvvz0NprbdG21BEdFLSHtF+vL2FpexzduLGBEUrzXcUQuyJCEWL5963T2+k/wkzUlXseRMKGilpC1+0g9P1hdzLKpWSyfNcrrOCJ9YnF+Oh+ePZqH15axraLO6zgSBlTUEpJa2jv4+99vJWVQHN++dbqumZaI8o0bC0gfksAXntyqs8DlnFTUEpJ++FIJRZUN/Odt0zXkLRFn6OA4fvDRmeyrPcG/P7/b6zgS4lTUEnI27T/Kz9ft5Y65OSyZnOl1HJF+MX9iGp9dOJ4nNh7kpV1VXseREKailpBS19TKF57cyujhg/iazvKWCPcPS/MpGJnCV57ZRnVDs9dxJESpqCVkOOf40v91fmD95I7ZDEkI6uZuImErIdbHf98+ixMt7fzDU+8SCDivI0kIUlFLyPjVhn28vLuKr143hVk5w7yOIzIg8jKT+cZNBawvqeGnr5Z6HUdCkIpaQsI7B4/xnReKWFqQyV8vGOd1HJEB9fG5Y7jlomx++HIxG0pqvI4jIUZFLZ6ra2rl/ifeIWtoIt+7baYuxZKoY2Z865Zp5GUM4fNPvsOR4ye9jiQhREUtnmrvCHD/E+/gb2jhpx+fzdDBcV5HEvHE4PhYHrrzYlraOrj/iXdo6wh4HUlChIpaPPUfLxSxobSGf79lmo5LS9SbmDGE73x4BpsPHONfVuzEOZ1cJqDTasUzT2+u4Fcb9vHp+eP46Jwcr+OIhISbZo5ix+Hj/HxtGfmZyXxq/jivI4nHtEctnnjn4DH+8dntzJ+QytdumOJ1HJGQ8uVrJ3P1lAz+7bldrC/xex1HPKailgFXfrSJe/53M5kpCTz48dnE+fTPUKQ7X4zxo9svIi9jCPf9dgt7/Y1eRxIP6RNSBlRdUyuffuxtWto6ePRTlzBc83iL9GhIQiy/+OQc4n0xfOZ/NlHT2OJ1JPGIiloGTHNbB599vJDyoyf5xSfnkJeZ7HUkkZCWM2Iwj3xyDpX1zdz1P5s40aI7bUUjFbUMiEDA8cWntrJp/zF+8LGZzBuf6nUkkbBw8djhPPjx2ew8XM+9v9lMa7su24o2Kmrpd845vrFiByu3V/K166dw44xRXkcSCStXTcnkP26dzvqSGv6/pzUneLTR5VnSr5xzfHvlbn7z1kH+ZvF47l6Y63UkkbD00Tk5+Bta+N6qPQwdFMe/3jxVs/hFCRW19KsfvlzCL9bv41OXjeWBZZP1wSJyAe67YgL1J9v4+boyYmNi+Kcbp+hnKgqoqKXf/Oy1vfz4lRI+cvFo/vkm/fYvcqHMjAeum0xrR4BHX99HnK/zsX62IpuKWvrFT14p4b9eKuammaP4zodnEBOjDxKRvmBmfOPGAto7XOeetc/40tJJKusIpqKWPuWc4/ur9/Dgq3u55aJsvnfbDHwqaZE+ZWb8681TaQ84Hnx1LydbAxoGj2Aqaukzzjn+/fnd/GrDPu6Ym8O3PjRde9Ii/SQmxvjWh6aRGBfDo6/v40RLO9++dbp+MY5AKmrpE20dAb7+7A5+X1jOp+eP459vKtBv9yL9LCamcxg8OTGOH79SQmNrOz/86CziY3XlbSRRUcsFO9HSzn2/3cLaYj+fXzKRv78mXyUtMkDMjC9ek09KYiz//vxujje18dAnZpOSqHu7R4qgfu0ys2VmtsfMSs3sgR5eNzP7cdfr28xsdrfX9pvZdjPbamaFfRlevFfd0MzHHnmTDaU1/Met0/miTmoR8cTdC8fz/Y/M5K2yWm772RtUHGvyOpL0kXMWtZn5gAeB64AC4A4zKzhjseuAvK4/9wA/O+P1K51zs5xzcy48soSK0upGbn3oDfZWn+AXn7yYO+aO8TqSSFS77eLRPH7XXI4cb+aWh95gW0Wd15GkDwSzRz0XKHXOlTnnWoEngeVnLLMceNx1egsYZmYj+zirhJCXd1XxoQdfp7mtgyfvuZQlkzO9jiQiwPyJafzhb+eTEBvDx37+Fi/uqPQ6klygYIo6Gyjv9rii67lgl3HAajPbbGb3nG9QCQ2BgOPHr5Rw9+OFjEsbzJ/uv5yZOcO8jiUi3eRlJvPsfQuYlJXMvb/ZzH++WESH5gcPW8EUdU8HHM/8G/+gZRY452bTOTz+OTNb1OObmN1jZoVmVuj3+4OIJQOtseuksR+8VMyHZo3i6Xvnkz1skNexRKQH6ckJPHnPpdwxN4eHXtvLpx97m2MnWr2OJechmKKuAHK6PR4NHA52Gefcqf9WA8/SOZT+Ps65R5xzc5xzc9LT04NLLwNm1+F6bv7pBlbvquTrN0zhhx+bRWKcz+tYIvIBEuN8/MetM/jOrdPZWHaUG3+yge0Vx72OJb0UTFFvAvLMLNfM4oHbgRVnLLMC+GTX2d+XAsedc0fMLMnMkgHMLAlYCuzow/zSz5xzPP7mfj700Os0Nrfzm7vncffC8TqzWySM3D53DE/dexnOOW792ev8Yl2ZbpUZRs55HbVzrt3M7gdWAT7gUefcTjO7t+v1h4GVwPVAKdAE/HXX6pnAs10f6rHAE865F/t8K6RfHG9q48vPvMuqnVVcMSmd739kJmlDEryOJSLnYVbOMJ7//EK+8sw2vrVyN+tK/PzXR2eSkZzodTQ5B3Mu9H6rmjNnjiss1CXXXnptTzUPPLOdmsYWvrJsMp+5PFfTgYpEAOccT7x9kG8+t4uk+Fi+++EZXF2gqza8Zmabz3YJs+aZk/doaG7jgWe28enHNpGcGMsf7pvPZxeNV0mLRAgz4855Y/nz/ZeTkZLI3Y8X8ve/36oTzUKYphCV09aX+Hngme0cOX6SexdP4AtX5+mEMZEIlZeZzJ8+t4CfvlrKQ6+Wsr6khn//0FSWTdMUGKFGe9RCVX0z9z+xhb/61dskxMXw9N/O54HrJqukRSJcfGwMX7wmnxX3X05mSgL3/mYLf/ubzRyuO+l1NOlGx6ijWHtHgMffPMAPXiqmtSPA566YyN8sHq+CFolCbR0BHllXxk/WlGAY/++qidx9+XjdiWuAfNAxahV1lNpYVsu/PbeLnYfrWZSfzr/dPJVxaUlexxIRj5UfbeKbz+1i9a4qxqcl8S83T2VRvua26G8qajltr7+R77xQxEu7qshKSeSfbizg+ulZui5aRN7j1T3V/OuKneyvbWLJ5Ay+smwyk7KSvY4VsVTUQm1jC//9Sgm/3XiQxNgY7rtyInctyGVQvIa5RaRnzW0dPPb6fh56rZQTLe18ePZovrg0n5FDNXVwX1NRR7GjJ1r55foyfv3GfprbA9wxN4e/uyqf9GRNXCIiwTl2opUHXy3l8TcPYAafnj+OexaNJ1UTIPUZFXUU6l7QTW0d3DB9JF+4Oo+JGRq6EpHzU360iR+8VMwftx4iMdbHJy4dw2cXjdfsZn1ARR1FquubeeyN/TzeraA/f1Ue+ZkqaBHpG6XVDTz46l7+tPUQcb4Y7pg7hr9ZPF5D4hdARR0Fiirr+eX6ffxp6yHaA47rp4/k71TQItKP9tec4KHXSvnDlkOYwU0zRnHX5blMyx7qdbSwo6KOUIGAY31pDb9cX8b6khoGxfn46JzR3HV5LmNTdamViAyM8qNN/GrDPv6vsJwTrR3MzR3BZy7P5eopmfg0/XBQVNQRpqaxhf8rrODJTQc5UNtEenICn54/jjvnjWHY4Hiv44lIlKpvbuOpTeU89vp+DtWdJGfEIG6/ZAwfmTNax7HPQUUdAZxzvFlWyxMbD7JqZyVtHY6540bw8XljuG56FgmxusxKREJDe0eA1buq+PUb+9m47yixMcZVUzK4fe4YFuWlay+7Bx9U1LopR4gr8zfyx62H+eM7hzh4tImUxFg+celYPj53DHk6/iwiISjWF8P100dy/fSRlPkb+f2mcp7eXMGqnVVkDxvErbOzWT4rm4kZQ7yOGha0Rx2Cahtb+PO7h3l262HeLa/DDBZMSOOWi7K5YcZIzcUtImGntT3Ay7ureHJTORtK/AQcTB2VwvJZo7hp5qioP2NcQ99hoKq+mdU7K3lxZyVvlR2lI+AoGJnCLRdlc9PMUWQN1fEdEYkM1fXNPLftCH969y87I/NyR3DdtJFcU5DJqGHRV9oq6hB1oPYEq3ZW8uKOSrYcrANgQnoSy6ZlcfPMbM2rKyIRb1/NCVZsPcyftx2mtLoRgBmjh7K0IJOlU7PIyxgSFfciUFGHiOa2Dt7ed5TX9vhZW1zNXv8JAKZnD2XZtCyunZqpmcNEJGrt9TeyemcVq3ZWsrW8DoDctCSumJTOovx0Ls1Njdj7E6ioPeKcY1/NCdYV+1lb7OfNslqa2wLEx8Zw6fhUrshP55qCTHJGDPY6qohISKmqb2b1ripe2lXFxrJaWto7PzvnjhvBovw0FuWnMykzOWL2tlXUA8Q5R1nNCTaWHeWtslreKquluqEF6PytcHF+OosnRfZvhSIifa25rYON+46yrtjPumI/JV1D5OnJCczLHcG88anMyx0R1sPkujyrn7R3BCiqbOCd8jre3tdZzv6uYs5ITuDS8alcOj6VBRNTNVOYiMh5Sozzde7o5KcDcLjuJOtL/LxeWsvGfbU8t+0IACOS4rlk3HDm5aYyN3cEk7KSifPFeBm9T2iPuhcqjzfzzsFjbC2v452DdWw/dJyTbR0AZKb8pZgvHZ/KuNTBYfubnYhIuHDOcfBoExv3HWVj2VE27qul4thJABLjYpg2aiizcoYxM2cYs3KGMXr4oJD8bNbQdy855zhUd5Jdh+vZdaSeXYfr2X7oOEeONwMQ74uhYFQKF43p/Iu/KGc4OSNC8y9fRCTaHKo7SeH+o7xbfpyt5cfYcbie1vYAAGlD4pk5ehjTRw9lysgUCkamhER5a+j7AzS3dbDX30jRkYbTpbzrSD3HT7YBYAbj05K4ZNyI08VcMCpFU3aKiISo7GGDyJ7VOfsZdE62sqeyga3lx9jaVd5r9lRzaj81OTGWKVkpFIxKYcrIZKaMTCE/MzlkJpeKmj3qhuY2SqsbT/8p6fpv+bGm039ZCbExTB6ZwtRRnb9lFYxKYXJWMoPjo/73GRGRiNLU2k5RZQO7u3bQdh+pp6iygabWzsOZZpAzfDATM4aQlzGECRlDmNj1JyUxrs/zRPUe9d2/LmTHoeNU1jeffi7eF8P49CRmjB7KrbOzyctIJj9zCLlpScRGwIkHIiLywQbHxzJ7zHBmjxl++rlAwHHgaBO7j9Szp7KBUn8jpVWNbCipobUjcHq5zJQEvrR0Eh+ZkzMgWSO+qJMSfMyfkMrEzCHkZSQzMWMIOcMHqZBFROQ9YmKM3LQkctOSuH76yNPPt3cEKD92sms0toHS6sYBndY5aoa+RUREQtUHDX1rt1JERCSEqahFRERCmIpaREQkhKmoRUREQpiKWkREJISpqEVEREJYUEVtZsvMbI+ZlZrZAz28bmb2467Xt5nZ7GDXFRERkbM7Z1GbmQ94ELgOKADuMLOCMxa7Dsjr+nMP8LNerCsiIiJnEcwe9Vyg1DlX5pxrBZ4Elp+xzHLgcdfpLWCYmY0Mcl0RERE5i2CKOhso7/a4ouu5YJYJZl0AzOweMys0s0K/3x9ELBERkcgXTFH3dJPOM+cdPdsywazb+aRzjzjn5jjn5qSnpwcRS0REJPIFc1OOCqD7LUJGA4eDXCY+iHVFRETkLILZo94E5JlZrpnFA7cDK85YZgXwya6zvy8FjjvnjgS5roiIiJzFOfeonXPtZnY/sArwAY8653aa2b1drz8MrASuB0qBJuCvP2jdc73n5s2ba8zswHluU0/SgJo+/H5eipRtiZTtAG1LqNK2hJ5I2Q7o+20Ze7YXQvI2l33NzArPdvuwcBMp2xIp2wHallClbQk9kbIdMLDbopnJREREQpiKWkREJIRFS1E/4nWAPhQp2xIp2wHallClbQk9kbIdMIDbEhXHqEVERMJVtOxRi4iIhKWoK2oz+5KZOTNL8zrL+TCzb3bdoWyrma02s1FeZzpfZvY9Myvq2p5nzWyY15nOl5l9xMx2mlnAzMLurNZIusudmT1qZtVmtsPrLBfCzHLM7FUz2931b+vvvM50vsws0czeNrN3u7blX73OdCHMzGdm75jZcwPxflFV1GaWA1wDHPQ6ywX4nnNuhnNuFvAc8A2P81yIl4BpzrkZQDHwVY/zXIgdwK3AOq+D9FYE3uXuf4BlXofoA+3APzjnpgCXAp8L47+XFmCJc24mMAtY1jU5Vrj6O2D3QL1ZVBU18EPgy5xlvvFw4Jyr7/YwifDeltXOufauh2/ROcVsWHLO7XbO7fE6x3mKqLvcOefWAUe9znGhnHNHnHNbur5uoLMYerypUajrurNiY9fDuK4/YfnZZWajgRuAXw7Ue0ZNUZvZzcAh59y7Xme5UGb2LTMrB+4kvPeou7sLeMHrEFEq6LvciTfMbBxwEbDR4yjnrWu4eCtQDbzknAvXbfkRnTt8gYF6w2BuyhE2zOxlIKuHl74G/COwdGATnZ8P2g7n3J+cc18DvmZmXwXuB/55QAP2wrm2pWuZr9E5zPfbgczWW8FsS5gK+i53MvDMbAjwDPCFM0bUwopzrgOY1XUuyrNmNs05F1bnEZjZjUC1c26zmV0xUO8bUUXtnLu6p+fNbDqQC7xrZtA5xLrFzOY65yoHMGJQzrYdPXgCeJ4QLupzbYuZfQq4EbjKhfi1gr34ewk3wdwhTzxgZnF0lvRvnXN/8DpPX3DO1ZnZa3SeRxBWRQ0sAG42s+uBRCDFzH7jnPtEf75pVAx9O+e2O+cynHPjnHPj6Pxgmh2KJX0uZpbX7eHNQJFXWS6UmS0DvgLc7Jxr8jpPFNNd7kKQde5V/ArY7Zz7gdd5LoSZpZ+6qsPMBgFXE4afXc65rzrnRnf1yO3Amv4uaYiSoo4w3zGzHWa2jc6h/LC9ZAP4KZAMvNR1udnDXgc6X2Z2i5lVAJcBz5vZKq8zBavrhL5Td7nbDTwVzF3uQpWZ/Q54E5hkZhVm9hmvM52nBcBfAUu6fj62du3JhaORwKtdn1ub6DxGPSCXNkUCzUwmIiISwrRHLSIiEsJU1CIiIiFMRS0iIhLCVNQiIiIhTEUtIiISwlTUIiIiIUxFLSIiEsJU1CIiIiHs/we9CmyHUkG+cwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t分布概率密度函数的值（x0=2，df=2）： 0.06804138174397717\n",
      "\n",
      "t分布在x0=1.5 处的累积分布函数值： 0.863803437554497\n",
      "\n",
      "使用Scipy计算概率密度函数的值： 0.06804138174397717\n",
      "\n",
      "使用Scipy计算累积分布函数的值： 0.8638034375544996\n"
     ]
    }
   ],
   "source": [
    "# 应用很广泛的t分布（学生t分布）\n",
    "\n",
    "import scipy.special as ss\n",
    "\n",
    "# t分布的概率密度函数\n",
    "# d 自由度\n",
    "def tf(x, d):\n",
    "    return ss.gamma((d+1)/2)*((1+(x**2)/d)**(-(d+1)/2))/((d*np.pi)**0.5*ss.gamma(d/2))\n",
    "\n",
    "d = 2\n",
    "x = np.arange(-4, 4, 0.01)\n",
    "t = ss.gamma((d+1)/2)*((1+(x**2)/d)**(-(d+1)/2))/((d*np.pi)**0.5*ss.gamma(d/2))\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(x, t)\n",
    "plt.show()\n",
    "\n",
    "# 手工计算的概率密度函数和分布函数的值\n",
    "# tf(随机变量, 自由度)\n",
    "print('t分布概率密度函数的值（x0=2，df=2）：', tf(2, 2))\n",
    "\n",
    "# 求概率密度函数的积分即可求出分布函数的值\n",
    "# integ.quad 求积分\n",
    "# integ.quad(分积分函数，积分下限，积分上限，args=自由度)\n",
    "print('\\nt分布在x0=1.5 处的累积分布函数值：', integ.quad(tf, float('-Inf'), 1.5, args=(2))[0])\n",
    "\n",
    "# F(x)满足 自由度为2的t分布\n",
    "\n",
    "# 通过Scipy 函数计算二者的值，其他函数的调用和前面正态分布类似，只不过t多了一个自由度参数\n",
    "# pdf(随机变量，自由度) => 计算概率密度函数\n",
    "# 与上面tf(2, 2) 计算结果对比\n",
    "print('\\n使用Scipy计算概率密度函数的值：', st.t.pdf(2, 2))\n",
    "# x0 = 1.5 自由度为2\n",
    "# 与上面integ.quad(tf, float('-Inf'), 1.5, args=(2))[0] 计算结果对比\n",
    "print('\\n使用Scipy计算累积分布函数的值：', st.t.cdf(1.5, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02633acb",
   "metadata": {},
   "source": [
    "#### 伽马分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cbfb2a1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.14936120510359183\n",
      "0.14936120510359185\n"
     ]
    }
   ],
   "source": [
    "# 伽马分布是一个连续分布，伽马分布的函数在Scipy.stats包\n",
    "# 伽马函数在scipy.special包\n",
    "# 示例：当x=3，a=2 时，伽马分布的概率密度函数值是多少？\n",
    "\n",
    "x, a = 3, 2\n",
    "gamma_pdf = x**(a-1)*np.exp(-3)/ss.gamma(2)\n",
    "\n",
    "# 比较手工计算与stats包gamma分布的pdf函数计算的值\n",
    "# 二者相等。后面一个值是调用pdf函数计算的\n",
    "print(gamma_pdf)\n",
    "print(st.gamma.pdf(3, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1d2b732",
   "metadata": {},
   "source": [
    "#### Beta分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "00983476",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4374000000000001\n",
      "0.4373999999999995\n"
     ]
    }
   ],
   "source": [
    "# Beta分布概率密度函数\n",
    "# a = 3, b = 4, x = 0.1\n",
    "\n",
    "beta_pdf = 0.1**2*0.9**3/ss.beta(3, 4)\n",
    "\n",
    "print(beta_pdf)\n",
    "print(st.beta.pdf(0.1, 3, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bc06521",
   "metadata": {},
   "source": [
    "#### Scipy各种概率分布对应表\n",
    "\n",
    "|  概率分布   | Scipy对应的名称 | 描述 |\n",
    "|  ----  | ----  | ----  | \n",
    "| beta  | beta | 贝塔分布 |\n",
    "| binomial  | binom | 二项分布 |\n",
    "| Cauchy  | cauchy | 柯西分布 |\n",
    "| chi-squared  | chi2 | 卡方分布 |\n",
    "| exponential  | expon | 指数分布 |\n",
    "| F  | f | f分布 |\n",
    "| gamma  | gamma | 伽马分布 |\n",
    "| geometric  | geom | 几何分布 |\n",
    "| hypergeometric  | hypergeom | 超几何分布 |\n",
    "| log-normal  | lognorm | 对数正态分布 |\n",
    "| logistic  | logistic | 逻辑斯蒂分布 |\n",
    "| negative binomial  | nbinom | 负二项分布 |\n",
    "| normal  | norm | 正态分布 |\n",
    "| Possion  | poisson | 泊松分布 |\n",
    "| Student's t  | t | （学生）t分布 |\n",
    "| uniform  | uniform | 统一分布 |\n",
    "| Weibull  | exponweib | 威布尔分布 |\n",
    "| Wilcoxon  | wilcoxon | 威尔科克孙分布 |\n",
    "\n",
    "\n",
    "> 调用方式\n",
    "```\n",
    "scipy.stats.norm.pdf()\n",
    "scipy.stats.norm.cdf()\n",
    "scipy.stats.norm.sf()\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4307f36e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
